// RSsimulation.cpp : Defines the entry point for the console application.
//
#include <time.h>
#include <math.h>
#include <stdlib.h>

#ifdef __cplusplus
extern "C" {  
#endif	

void rsdepth(double *points, double *query,double *depth, int*sz);

#ifdef __cplusplus
}
#endif

typedef struct __MyPoint
{
	double x,y,angle;
//	int operator ==(MyPoint two)
//	{
//		if(two.x==x &&two.y==y)
//			return 1;
//		return 0;	
//	}	
}MyPoint;

double minimum(double a,double b)
{
	if (a<b) return a;
	return b;
}

double maxmum(double a,double b)
{
	if (a>b) return a;
	return b;
}

int compare (const void * a, const void * b)
{
  return ( ((MyPoint*)a)->angle - ((MyPoint*)b)->angle );
}

int cmp (const void * a, const void * b)
{
  return ( *(double*)a - *(double*)b );
}

typedef struct __MyLine
{
	MyPoint a,b;
}MyLine;

int numberOfPoints = 9;
MyPoint q,*P;

double Depth;
int algo;
int did2,did1;
int m=1,n=9;
char ch;

double RSDepth(double qx, double qy);
double _RSDepth(MyPoint q);

//.C("rsdepth", points=as.double(z), query=as.double(q),  depth=as.double(dp))

void rsdepth(double *points, double *query,double *depth, int *sz)
{	
	int i;
	int count=0;
	P = (MyPoint*) (malloc(sizeof(MyPoint)*300));
	if (P==NULL) exit (1);
	n=sz[0];
	numberOfPoints=n;
	for( i=0;i<numberOfPoints;i++)	{
		P[i].x=points[count++];
		P[i].y=points[count++];
	}	
/*
	for( i=0;i<numberOfPoints;i++)	{
		P[i].y=points[count++];
	}	
*/
	depth[0]= RSDepth(query[0],query[1]);
	return;
}

double _RSDepth(MyPoint q)
{
	double angle[300];
	
	int total[300];
	
	int endOf[300];
	int i;
	int count = 0;//points in between 
    int index=0;//of points for which we have found antipodals
    double minDepth = 100000;    
	double ray[100];
	double current_angle=0,largest_angle = 0;
    
	for(i=0;i<numberOfPoints;i++)
	{
		double slope ;
		if(q.y-P[i].y==0&&(q.x-P[i].x)==0)
			slope=0;
		else if ((q.x-P[i].x)==0)
			slope = 99999999999999999999.0;
		else
		slope = (q.y-P[i].y)/(q.x-P[i].x);
		angle[i]=atan(slope)*57.35;

		if (angle[i] < 0)   {
            angle[i] += 180;                
        }

        if (q.y > P[i].y) {
            angle[i] += 180;                 
        }	
		P[i].angle=angle[i];
	}

	//sort
	qsort( angle, numberOfPoints, sizeof(angle[0]), cmp);
	//sort(angle,angle+numberOfPoints);
	qsort( P, numberOfPoints, sizeof(P[0]), compare);
	//sort(P,P+numberOfPoints);


	//Fill array from number of points between a point and its antipodal
	
    
	for ( i = 0; index < numberOfPoints  ; i++)
    {
		count++;
        if ( (i-index)>= numberOfPoints ||angle[i % numberOfPoints ] > (angle[index] + 180) ||
			
			((angle[i % numberOfPoints ] < (angle[index]) && (angle[i % numberOfPoints ] + 360 > (angle[index] + 180)))))
        
		{
			count--;
        
			total[index] = count - 1;
            endOf[index] = i - 1;
            index++;
            count--;
            i--;

        }
    }

	//calculate depth for first ray between 1 and n points.
	
    ray[0] = 0;
	did2=0;
    for (i = 0; i <= total[0]; i++)
		ray[0] += numberOfPoints - ( (total[i])+ i +1);
    minDepth = ray[0];

	if(ray[0]>=1.0/9.0*numberOfPoints*numberOfPoints)
	{
		current_angle+=angle[0]+360-angle[numberOfPoints-1];
	}
	else;
		

	for( i=1;i<numberOfPoints+1;i++)
    {
        ray[i] = ray[i - 1] - numberOfPoints + 2*total[i-1]+1;

		if(ray[i]>=1.0/9.0*numberOfPoints*numberOfPoints)
		{
			current_angle+=angle[i]-angle[i-1];
		}
		else
		{
			if(current_angle > largest_angle)
			{
				largest_angle=current_angle;
				current_angle=0;			
			
			}
		}


		if (minDepth > ray[i])
		{
            minDepth = ray[i];            
			did2=i;
		}
    }
    return minDepth;

}

double RSDepth(double qx, double qy)
{
	MyPoint pt;
	pt.x=qx;
	pt.y=qy;
	return _RSDepth(pt);
}
int InBoundingBox(MyLine a, MyPoint p)
{
	
            
    if (p.x <minimum(a.a.x, a.b.x))
        return 0;

    if (p.x > maxmum(a.a.x, a.b.x))
        return 0;

    if (p.y <minimum(a.a.y, a.b.y))
        return 0;

    if (p.y > maxmum(a.a.y, a.b.y))
        return 0;

    return 1;      
        
}
MyPoint Intersects( MyLine first, MyLine second)
{
	
	MyPoint intersect;
	double a1 = (first.b.y - first.a.y);
	    
	double b1 = (first.a.x - first.b.x);
	
	double c1 = (first.b.x*first.a.y - first.a.x*first.b.y);
        	

	//Line to test against
	double a2 = (second.b.y - second.a.y);
	double b2 = (second.a.x - second.b.x);
	double c2 = (second.b.x*second.a.y - second.a.x*second.b.y);

	double denom = a1*b2 - a2*b1;

	intersect.x=intersect.y=0;
	
	//this line
                           
        	
	//Check for parallel lines
	if(denom == 0) { return intersect; }
        	
	//Get the intersection point
	intersect.x = ( (b1*c2 - b2*c1)/denom);
	intersect.y = ( (a2*c1 - a1*c2)/denom);

	//InBoundingBox tests to see if the point is in the bounding box for
		//a line segment. The point must be in both lines' bounding boxes to
		//register and intersection.
	if( InBoundingBox( first,intersect )==0 || InBoundingBox( second,intersect )==0 )
		{
			intersect.x=0;
			intersect.y = 0;
		}
            	
	return intersect;

}

int IsUnique(MyPoint p)
{
	int i;
    for(i=0;i<numberOfPoints;i++)    {
        if (p.x == P[i].x && p.y == P[i].y)
            return 0;
    }

    return 1;        
}


 MyPoint RSMedian()
{
	    double maxmum = 0;
    double t = 0;
    int i;        

	MyPoint dum,ok;
    MyLine l1;
	MyLine l2;

	l1.a=P[0];
	l1.b=P[1];    
	l1.a=P[2];
	l1.b=P[3];
    
    for (i = 0; i < numberOfPoints; i++)
    {
        int j;
		l1.a = P[i];
        for (j = i + 1; j < numberOfPoints; j++)
        {
			int k;       
            l1.b = P[j];
             
            for ( k = 0; k < numberOfPoints; k++)
            {
                 int l;    
                if (i == k || j == k)
                    continue;
                l2.a = P[k];
                for (l = k + 1; l < numberOfPoints; l++)
                {

                    if (i == l || j == l)
                        continue;
                    l2.b = P[l];
                    dum = Intersects(l1, l2);

					if (dum.x != 0 && dum.y != 0)  {                        
                        t =RSDepth(dum.x , dum.y);
                        if (t > maxmum && IsUnique(dum)) { maxmum = t;ok = dum;}
                    }

                }                    
            }
                
        }
            
    }
    
	return ok;        
}

